allocatemem; \\ doubles the memory allocation
allocatemem; \\ doubles the memory allocation
allocatemem; \\ doubles the memory allocation
allocatemem; \\ doubles the memory allocation
\p 100;
default(format, "f0.9");
defaultFormat = "f0.9";
defaultMinX = 1.0*10^(-50); 
defaultMaxX = 1.0*10^50;
defaultMinY = 1.0*10^(-22);
defaultMaxY = 1.0*10^(22);
maxX = defaultMaxX;
minX = defaultMinX;
minY = 1.0*10^(-10000); \\ 5
true = 1;
false = 0;
pdfs = 1;
directory = "Need to dynamically assign.";

lBracket1 = 0;
rBracket1 = 0;
lBracket2 = 0;
rBracket2 = 0; \\ THIS WAS rBracket3
fString1 = "";
fString2 = "";

\\ initialize pdfs
startMinY = 1.0*10^(-1000);
midMinY = 1.0*10^(-50);
endMinY = 1.0*10^(-22);
defaultMaxSearchX = 100000000000;

\\ Assigns strings and finds valid regions for each pdf. There are three
\\ searches for valid region. Each search uses the previous searches results
\\ to narrow down where most of the area is. This is required because if
\\ the region is allowed to be too big, then most of the integral is extreemly flat,
\\ and the integrator performs poorly. 
initialize(pdfString,pdfID) = { \\ ID is 1 or 2
	maxSearchX = defaultMaxSearchX;
	minY = startMinY; 
	lBracketTemp = getBoundary(pdfString,minX,1);
	rBracketTemp = getBoundary(pdfString,lBracketTemp,0);
	maxSearchX = min(rBracketTemp,maxSearchX);
	print("search1 done"," search distance: ",abs(lBracketTemp-rBracketTemp));
	print("lBracketTemp: ",lBracketTemp);
	print("rBracketTemp: ",rBracketTemp);
	
	minY = midMinY;
	lBracketTemp = getBoundary(pdfString,lBracketTemp,1);
	rBracketTemp = getBoundary(pdfString,lBracketTemp,0);
	maxSearchX = min(rBracketTemp,maxSearchX);
	print("search2 done"," search distance: ",abs(lBracketTemp-rBracketTemp));
	print("lBracketTemp: ",lBracketTemp);
	print("rBracketTemp: ",rBracketTemp);
	
	minY = endMinY;
	lBracketTemp = getBoundary(pdfString,lBracketTemp,1);
	rBracketTemp = getBoundary(pdfString,lBracketTemp,0);
	maxSearchX = min(rBracketTemp,maxSearchX);
	print("search3 done"," search distance: ",abs(lBracketTemp-rBracketTemp));
	print("lBracketTemp: ",lBracketTemp);
	print("rBracketTemp: ",rBracketTemp);
	
	if(pdfID == 1,
		fString1 = pdfString;
		lBracket1 = lBracketTemp;
		rBracket1 = rBracketTemp
		,
		fString2 = pdfString;
		lBracket2 = lBracketTemp;
		rBracket2 = rBracketTemp
	);
}

\\ integrate from a to b, for one or two pdfs
integrate(a,b,withX) = {
	newLeft = max(a,lBracket1);
	newRight = min(b,rBracket1);
	aPhrase = getPhrase(fString1, newLeft,newRight,1);
	bPhrase = getPhrase(fString1,newLeft,newRight,0);
	if(withX,
		area1 = intnum(x = aPhrase, bPhrase, evalFunction(fString1,x)*x);
		,
		area1 = intnum(x = aPhrase, bPhrase, evalFunction(fString1,x));
		
	);
	if(pdfs==1,
		return(area1);
	);
	newLeft = max(a,lBracket2);
	newRight = min(b,rBracket2);
	aPhrase = getPhrase(fString2, newLeft,newRight,1);
	bPhrase = getPhrase(fString2,newLeft,newRight,0);
	if(withX,
		area2 = intnum(x = aPhrase, bPhrase, evalFunction(fString2,x)*x);
		,
		area2 = intnum(x = aPhrase, bPhrase, evalFunction(fString2,x));
	);
	return(area1 + area2);
}

\\ helper function used in several places in the integrate method
bPhrase2 = 0;
getPhrase(fString5, newLeft,newRight,isA) = {
	temp = 1000000;
	delta = 1/temp;
	outLeft = evalFunction(fString5,newLeft);
	outLeft2 = evalFunction(fString5,newLeft+delta);
	slope = (outLeft2-outLeft)/delta;
	if(isA,
		return(if(slope < -10000,[newLeft,-.9],newLeft)); \\ different parameters if extreem at left
		,
		return([newRight,bPhrase2]);
		\\ return(if(newRight > 10^10,[newRight,-.1],newRight));
		\\ return(if(newRight > 10^10,[+1],newRight)); \\ integrates better to infinity in some cases
	);
}

calculateAllVARCTEs() = {	
	foundArea = integrate(minX,maxX,0);
	foundAreaInBounds = integrate();
	default(format, "f0.16");
	foundAreaString = Str(foundArea);
	if(foundArea < .0000000000000001,
		foundAreaString = "0";
	);
	default(format, "f0.10");
	print1(fString1);
	writeToFile("<FUNCTION_STRING>");
	writeToFile(fString1);
	print("minX: ",minX);
	print("maxX: ",maxX);
	if(pdfs==2,
		print(Str(" + ", "  TOTAL AREA: ",foundAreaString));
		print(fString2);
		writeToFile(fString2);
		writeToFile("</FUNCTION_STRING>");
		print("bounds: ", "    ", lBracket1,"     ",rBracket1);
		print("bounds: ", "    ", lBracket2,"     ",rBracket2);
		,
		print(Str("  TOTAL AREA: ",foundAreaString));
		print("bounds: ", "    ", lBracket1,"     ",rBracket1);
		writeToFile("</FUNCTION_STRING>");
	);
	writeToFile("");
	writeToFile("<TOTAL_AREA>");
	writeToFile(foundAreaString);
	writeToFile("</TOTAL_AREA>");
	writeToFile("");
	writeToFile("<BOUNDS>");
	writeToFile(Str(lBracket1,"     ",rBracket1));
	if(pdfs==2,
		writeToFile(Str(lBracket2,"     ",rBracket2));
	);
	writeToFile("</BOUNDS>");
	writeToFile("<AREA_LEFTCTE_VAR_RIGHTCTE>");
	resetGetVAR2();
	print(foundAreaString);
	
	print("AREA AT VAR        LEFT CTE         VAR          RIGHT CTE");
	calculateVARCTE(.01);
	calculateVARCTE(.05);

	calculateVARCTE(.10);
	calculateVARCTE(.15);
	calculateVARCTE(.20);
	calculateVARCTE(.25);
	calculateVARCTE(.30);
	calculateVARCTE(.35);
	calculateVARCTE(.40);
	calculateVARCTE(.45);
	
	calculateVARCTE(.50);
	
	calculateVARCTE(.55);
	calculateVARCTE(.60);
	calculateVARCTE(.65);
	calculateVARCTE(.70);
	calculateVARCTE(.75);
	calculateVARCTE(.80);
	calculateVARCTE(.85);
	calculateVARCTE(.90);
	
	calculateVARCTE(.95);
	calculateVARCTE(.99);
	\\ print("");
	
	writeToFile("</AREA_LEFTCTE_VAR_RIGHTCTE>");
	resetBounds();
}

\\ helper function
evalFunction(fString,x) = {
	xVar = x;
	out = eval(fString);
	return(out);
}

binarySearch_findEdge(fString,lftBrckt,rghtBrckt,leftEdge) = {
	done2 = 0;
	while(!done2,
		midpoint2 = (lftBrckt + rghtBrckt)/2;
		outputMid = evalFunction(fString,midpoint2);
		\\print("====brackets: ",lftBrckt,"  ,  ",rghtBrckt);
		
		if(abs(lftBrckt-rghtBrckt) < defaultMinX||lftBrckt==midpoint2||rghtBrckt==midpoint2,
			done2=1;
		);
		
		if(outputMid > minY, \\ use default to prevent excessive search
			if(leftEdge==true,
				
				rghtBrckt = midpoint2;,
				lftBrckt = midpoint2;
			);
			,
			if(leftEdge==true,
				lftBrckt = midpoint2;,
				rghtBrckt = midpoint2;
			);
		);
	);
	
	\\ return slightly inside boundaries to ensure positive function
	if(leftEdge==true,
		print("BOUNDARY FOUND, EVALUATION: ",evalFunction(fString,rghtBrckt));
		return(rghtBrckt);
		,
		print("BOUNDARY FOUND, EVALUATION: ",evalFunction(fString,lftBrckt));
		return(lftBrckt);
	);
}

defaultScaleFactor = .1;
scaleFactor = defaultScaleFactor;
minScaleFactor = scaleFactor * .1;
multiplier = 0;
leftFound = 0;
startWidth = .0001;

\\ this finds the left and right boundaries of the pdf (where the output is above a certain threshold). This is
\\ necessary because peri integration is not dynamic, and therefore performs better when this 
\\ is specified exactly. 
getBoundary(fString,startSearch,leftBoundary) = {
	\\ entireWidth = abs(startSearch-max(rBracket1,rBracket2));
	width_2 = startWidth;
	\\ width_2 = abs(startSearch-max(rBracket1,rBracket2)) * .0001; \\ start a small proportion of entire width each time
	scaleFactor = defaultScaleFactor;
	multiplier = 1 + scaleFactor;
	if(leftBoundary,
		print(fString);
	);
	done_2 = 0;
	if(leftBoundary,leftFound = 0;,);
	leftBracket_2 = startSearch;
	rightBracket_2 = startSearch+width_2;
	count4 = 0;
	print("getBoundary start: leftBracket_2: ",leftBracket_2,"  rightBracket_2: ",rightBracket_2);
	
	while(!done_2,
		count4++;
		output_adjstBnds = evalFunction(fString,rightBracket_2);	
		endSearch = if(leftBoundary,output_adjstBnds > minY,output_adjstBnds < minY);
		if(endSearch,
			done_2 = 1;
			if(leftBoundary,
				leftFound = 1;
			);
			outputLeft =  evalFunction(fString,leftBracket_2);
			result4 = binarySearch_findEdge(fString,leftBracket_2,rightBracket_2,leftBoundary);
			print("result of binary search: ",result4, " leftEdge: ", leftBoundary);
			return(result4);
		);
		
		if(leftBoundary,
			if(leftBracket_2 > maxSearchX,
				scaleFactor *= .1; \\ if failed to find any valid region, make more precise
				if(scaleFactor < minScaleFactor,
					return(startSearch);
				);
				multiplier = 1 + scaleFactor;
				width_2 = startWidth;
				leftBracket_2 = startSearch;
				rightBracket_2 = leftBracket_2+width_2;
				print(" reducing width"," search distance: ",abs(startSearch-maxSearchX)," brackets: ","leftBracket_2:",leftBracket_2,"rightBracket_2:",rightBracket_2);
				next();
			);
			,
			if(leftBracket_2 > defaultMaxX,
				return(defaultMaxX);
			);
		);
		
		if(!leftBoundary,
			width_2 *= 1.1;
		);
		width_2 *= multiplier;
		leftBracket_2 = rightBracket_2;
		rightBracket_2 = startSearch + width_2;
		if(count4%1 == 0,
			\\ print("scaleFactor:",scaleFactor,"  leftBracket_2: ",leftBracket_2,"  leftBoundary: ", leftBoundary,"  maxSearchX",maxSearchX, "   ",fString);
		);
	);
}

\\ helper function
resetBounds() = {
	minX = defaultMinX;
	maxX = defaultMaxX;
}

\\ helper functions: some exponential and gamma functions are too large for perry,
\\ therefore, this prevents exceptions and wasted time
validExp(num) = {
	tempvalidexp = abs(num);
	if(tempvalidexp > 100000000,
		return(0);
	);
	return(1);
}
validGamma(num) = {
	if(num <= 0,
		return(0);
	);
	if(num > 100000,
		return(0);
	);
	return(1);
}
fReplace(x1) = {
	return(0);
}

\\ NOT USED.
buildStringMixed(fstring1,fstring2,args1,args2,scalar1,scalar2) = {
	fstring1 = Str(scalar1,"*",fstring1,"(xVar,",args1,")");
	fstring2 = Str(scalar2,"*",fstring2,"(xVar,",args2,")");
	mixedString = Str(fstring1,"+",fstring2);
	return(mixedString);
}

\\ NOT USED
buildString(fstring1,args1,scalar) = {
	fstring1 = Str(scalar,"*",fstring1,"(xVar,",args1,")");
	return(fstring1);
}

numCalculations = 0;
writeToFile(fileString) = {
	outputName = Str(directory,"/","PERI_OUTPUT.txt");
	write(outputName,fileString);
}

\\ writes to a text file that is intermittently checked by AMOOF to indicate when
\\ peri is done or failed
writeScriptDone(fileName4) = {
	write(fileName4,"done");
}
writeScriptFailed(fileName4) = {
	write(fileName4,"failed");
}

\\ find a single VARCTE according to a specific area
calculateVARCTE(desiredArea) = {
	vaR = getVar(desiredArea);
	foundArea = integrate(minX,vaR,0);
	cteLeft = getCte(vaR,1);
	cteRight = getCte(vaR,0);
	strng = Str(foundArea, " ; ", cteLeft, " ; ", vaR, " ; ", cteRight);
	writeToFile(strng);
	print(strng);
}

getCte(vaR,leftTail) = {
	if(leftTail,
			r3 = integrate(minX,vaR,1);
			r4 = integrate(minX,vaR,0);
			return(r3/r4 * 1.0);
			,
			r3 = integrate(vaR,maxX,1);
			r4 = integrate(vaR,maxX,0);
			return(r3/r4);
	);
}

\\ used to store previous search, thereby reduce 
\\ subsequent search space and improve speed
lastXFound = minX; 
resetGetVAR2() = {
	lastXFound = minX;
}

finalLeft() = {
	if(pdfs == 1,
			return(lBracket1);
			,
			return(MIN(lBracket1,lBracket2));
	);
}
finalRight() = {
	mx = max(1,4);
	if(pdfs == 1,
			return(rBracket1);
			,
			return(max(rBracket1,rBracket2));
	);
}

getVar(desiredArea) = {
	\\ leftBracket = max(finalLeft(),lastXFound);
	leftBracket = finalLeft();
	width = (finalRight()-leftBracket)/1000; 
	rightBracket = leftBracket + width;
	targetError = .0000001;
	done = 0;
	lastArea = 0;
	
	print("");
	print("expanding bounds");
	print("minX: ",minX);
	print("rightBracket: ",rightBracket);
	print("lBracket1: ",lBracket1);
	print("rBracket1: ",rBracket1);
	print("lBracket2: ",lBracket2);
	print("rBracket2: ",rBracket2);
	print("finalLeft: ",finalLeft());
	print("finalRight: ",finalRight());
	print("totalAreaInBounds: ",integrate(finalLeft(),finalRight(),0));
	print("bracketArea: ",integrate(leftBracket,rightBracket,0));
	print("lastArea: ",lastArea);
	
	print("--------expanding bounds");
	while(!done, 
		\\ print("-----in while loop");
		thisArea = integrate(leftBracket,rightBracket,0);
		print("thisArea: ", thisArea);
		if(thisArea <= lastArea && thisArea < desiredArea,
			print("------------EXPANSION NO LONGER IMPROVES");
			print("thisArea: ",thisArea);
			print("lastArea: ",lastArea);
			return(rightBracket); \\ expansion no longer improves, search has failed
			,
		);
		lastArea = thisArea;
		
		if(thisArea > desiredArea,
			done = 1;
			print("thisArea > desiredArea");
			,
			done = 0;
			width = width * 2;
			rightBracket = leftBracket + width;
		);
	);
	\\ narrow down
	done = 0;
	\\ print("bounds distance: ", abs(leftBracket-rightBracket));
	
	print("");
	print("-------------narrow down");
	baseLeft = leftBracket;
	while(!done,
		\\ print("-----in while loop");
		midpoint = (leftBracket + rightBracket)/2;
		thisArea = integrate(baseLeft,midpoint,0);
		print("thisArea: ",thisArea);
		areaDistance = abs(thisArea - desiredArea);
		if(thisArea > desiredArea,
			rightBracket = midpoint;,
			leftBracket = midpoint;
		);
		bracketWidth = abs(leftBracket - rightBracket);
		bracketWidthError = bracketWidth / midpoint * 1.0;
		areaError = areaDistance / 1.0;
		if(areaError < targetError,
			lastXFound = midpoint;
			return(midpoint * 1.0);
		);
	);
}

\\ NOT USED CURRENTLY. Uses interpolation to more quickly narrow down 
\\ the search space. Needs to be revised: sometimes hangs, also
\\ if the estimation is slightly on the side of the smaller half, then
\\ there is not much benefit for the extra calls required.
estimateMidpoint(leftBracket, rightBracket, desiredArea) = {
	lftB = leftBracket;
	leftArea = integrate(minX,lftB,0);
	rightArea = integrate(minX,rightBracket,0);		
	totalDifferenceArea = abs(leftArea - rightArea);
	desiredDifference = abs(leftArea - desiredArea);
	areaDiffProportion = desiredDifference / totalDifferenceArea;
	totalXDistance = abs(lftB - rightBracket);
	newDistance = areaDiffProportion * totalXDistance * 1.0;
	newLeftDist = newDistance;
	newRightDist = (totalXDistance - newLeftDist) * 1.0;
	finalX = lftB + newLeftDist;
	return(finalX * 1.0);
}

\\ START OF PDF FUNCTIONS
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\\ BETA ===================requires condition on x; difficulty finding area
Beta_(x,ALPHA,BETA,THETA) = { \\ single
	trap(,return(0),
		if(x<=0,return(0),);
		if(x>=THETA,return(0),);
		u1 = (x/THETA);
		if((1-u1)==0 && (BETA-1)<0,return(0),); \\ special case where results in division by zero
		num = gamma(ALPHA+BETA)*(u1*1.0)^ALPHA*(1.0-u1)^(BETA-1);
		denom = gamma(ALPHA)*gamma(BETA)*x;
		return(num/denom*1.0);
	);
}

\\ BURR
Burr_(x,ALPHA,GAMMA_1,THETA) = { \\ single
	trap(,return(0),
		num = ALPHA*GAMMA_1*(x/THETA)^(GAMMA_1);
		denom = x*(1.0+(x/THETA)^GAMMA_1)^(ALPHA+1);
		return(num/denom*1.0);
	);
}

\\ EXPONENTIAL
Exponential_(x,THETA) = { \\ single
	trap(,return(0),
		callcount++;
		if(!validExp(-x/THETA),
			return(fReplace(x));
		);
		num = exp(-x/THETA*1.0);
		denom = THETA;
		return(num/denom * 1.0);
	);
}

\\ GAMMA
Gamma_(x,ALPHA,THETA) = { \\ single
	ALPHA *= 1.0;
	THETA *= 1.0;
	trap(,return(0),
		if(!validExp(-x/THETA),
			return(fReplace(x));
		);
		if(!validGamma(ALPHA),
			return(fReplace(x));
		);
		num = (x/THETA*1.0)^ALPHA * exp(-x/THETA);
		denom = x * gamma(ALPHA);
		return(num/denom * 1.0);
	);
}

\\ GENERALIZED BETA ==============requires conditions on x; difficulty finding area
GeneralizedBeta_(x,ALPHA,BETA,TAU,THETA) = { \\ single
	trap(,return(0),
		if(x<=0,return(0),);
		if(x>=THETA,return(0),);
		u1 = (x/THETA)^TAU;
		if((1-u1)==0 && (BETA-1)<0,return(0),); \\ special case where results in division by zero
		num = gamma(ALPHA+BETA)*(u1)^ALPHA*(1-u1)^(BETA-1)*TAU;
		denom = gamma(ALPHA) * gamma(BETA)*x;
		return(num/denom*1.0);
	);
}

\\ GENERALIZED PARETO
GeneralizedPareto_(x,GAMMA_1,THETA) = { \\ single
	trap(,return(0),
		num = (1/GAMMA_1)*(THETA/GAMMA_1)^(1/GAMMA_1);
		denom = (x+THETA/GAMMA_1)^(1/GAMMA_1+1);
		return(num/denom);
	);
}

\\ INVERSE BURR
InverseBurr_(x,GAMMA_1,TAU,THETA) = { \\ single
	trap(,return(0),
		num = TAU*GAMMA_1*(x/THETA*1.0)^(TAU*GAMMA_1);
		denom = x * (1.0 + (x/THETA*1.0)^GAMMA_1)^(TAU + 1);
		return(num/denom);
	);
}

\\ INVERSE EXPONENTIAL
InverseExponential_(x,THETA) = { \\ single
	trap(,return(0),
		if(!validExp(-THETA/x),
			return(fReplace(x));
		);
		num = THETA*exp(-THETA/x);
		denom = (x*1.0)^2;
		return(num/denom * 1.0);
	);
}

\\ INVERSE GAMMA
InverseGamma_(x,ALPHA,THETA) = { \\ single
	trap(,return(0),
		if(!validExp(-THETA/x),
			return(fReplace(x));
		);
		num = (THETA/x)^ALPHA*exp(-THETA/x);
		denom = x*gamma(ALPHA);
		return(num/denom * 1.0);
	);
}

\\ INVERSE GAUSSIAN
InverseGaussian_(x,MU,THETA) = { \\ single
	trap(,return(0),
		z = (x-MU)/MU;
		if(!validExp((-THETA*(z*1.0)^2)/(2*x)),
			return(fReplace(x));
		);
		result = (THETA/(2*Pi*x^3))^(1/2)*exp((-THETA*z^2)/(2*x));
		return(result);
	);
}

\\ INVERSE PARALOGISTIC
InverseParalogistic_(x,TAU,THETA) = { \\ single
	trap(,return(0),
		num = (TAU*1.0)^2 * (x/THETA)^(TAU^2.0);
		denom = x * (1.0 + (x/THETA)^TAU)^(TAU + 1);
		return(num/denom * 1.0);
	);
}

\\ INVERSE PARETO
InversePareto_(x,TAU,THETA) = { \\ single
	trap(,return(0),
		num = TAU * THETA * x^(TAU - 1);
		denom = (x + THETA)^(TAU + 1);
		return(num/denom * 1.0);
	);
}

\\ INVERSE TRANSFORMED GAMMA
InverseTransformedGamma_(x,ALPHA,TAU,THETA) = { \\ single
	trap(,return(0),
		u1 = (THETA/x*1.0)^TAU;
		if(!validExp(-u1),
			return(fReplace(x));
		);
		num = TAU*u1^ALPHA*exp(-u1);
		denom = x*gamma(ALPHA);
		return(num/denom);
	);
}

\\ INVERSE WEIBULL =================not working: diverges
InverseWeibull_(x,TAU,THETA) = { \\ single
	trap(,return(0),
		if(!validExp(-(THETA/x)^TAU),
			return(fReplace(x));
		);
		num = TAU*(THETA/x*1.0)^TAU*exp(-(THETA/x*1.0)^TAU);
		denom = x;
		return(num/denom * 1.0);
	);
}

\\ LOGLOGISTIC 
Loglogistic_(x,GAMMA_1,THETA) = { \\ single
	trap(,return(0),
		num = GAMMA_1*(x/THETA)^GAMMA_1;
		denom = x*(1.0+(x/THETA)^GAMMA_1)^2;
		return(num/denom * 1.0);
	);
}

\\ LOGNORMAL ===============is this log or ln
logStor = sqrt(2*Pi);
Lognormal_(x,MU,SIGMA) = { \\ single
	trap(,return(0),
		z = (log(x)-MU)/SIGMA;
		if(!validExp(-(z^2/2)),
			return(fReplace(x));
		);	
		num = exp(-(z^2/2));
		denom = x * SIGMA * logStor;
		return(num/denom * 1.0);
	);
}

\\ PARALOGISTIC
Paralogistic_(x,ALPHA,THETA) = { \\ single
	trap(,return(0),
		num = ALPHA^2*(x/THETA)^ALPHA;
		denom = x * (1.0+(x/THETA)^ALPHA)^(ALPHA+1);
		return(num/denom * 1.0);
	);
}

\\ PARETO
Pareto_(x,ALPHA,THETA) = { \\ single
	ALPHA = ALPHA * 1.0;
	THETA = THETA * 1.0;
	trap(,return(0),
		num = ALPHA * THETA^ALPHA;
		denom = ((x + THETA)*1.0)^(ALPHA + 1.0);
		return(num/denom * 1.0);
	);
}

\\ SINGLE PARAMETER PARETO 
SingleParameterPareto_(x,ALPHA,THETA) = { \\ single
	trap(,return(0),
		if(x<=THETA,return(0),);
		num=ALPHA*(THETA*1.0)^ALPHA;
		denom = (x*1.0)^(ALPHA+1);
		return(num/denom);
	);
}

\\ TRANSFORMED BETA
TransformedBeta_(x,ALPHA,GAMMA_1,TAU,THETA) = { \\ single
	trap(,return(0),
		if(!validGamma(ALPHA+TAU),
			return(fReplace(x));
		);
		num = gamma(ALPHA + TAU) * GAMMA_1 * (x/THETA*1.0)^(GAMMA_1 * TAU);
		denom = gamma(ALPHA) * gamma(TAU) * x * (1.0 + (x/THETA)^GAMMA_1)^(ALPHA + TAU);
		return(num/denom);
	);
}

\\ TRANSFORMED GAMMA
TransformedGamma_(x,ALPHA,TAU,THETA) = { \\ single
	ALPHA *= 1.0;
	THETA *= 1.0;
	TAU *= 1.0;
	trap(,return(0),
		u1 = (x/THETA*1.0)^TAU;
		if(!validExp(-u1),
			return(fReplace(x));
		);
		num = TAU*(u1*1.0)^ALPHA*exp(-u1);
		denom = x*gamma(ALPHA);
		return(num/denom);
	);
}

\\ WEIBULL 
Weibull_(x,TAU,THETA) = { \\ single
	THETA *= 1.0;
	TAU *= 1.0;
	trap(,return(0),
		if(!validExp(-(x/THETA)^TAU),
			return(fReplace(x));
		);
		num = TAU*(x/THETA * 1.0)^(TAU)*exp(-(x/THETA* 1.0)^TAU);
		denom = x;
		return(num/denom* 1.0);
	);
}

print("");
print("");
/*
fstring = "Lognormal_(xVar,12.3254344,16.00231939)"
initialize(fstring,1);
foundArea = integrate(minX,maxX,0);
print(foundArea);
print("");

bPhrase2 = -.1;
foundArea = integrate(minX,maxX,0);
print("bPhrase2: ",bPhrase2);
print(foundArea);
print("");

bPhrase2 = -.01;
foundArea = integrate(minX,maxX,0);
print("bPhrase2: ",bPhrase2);
print(foundArea);
print("");

bPhrase2 = .1;
foundArea = integrate(minX,maxX,0);
print("bPhrase2: ",bPhrase2);
print(foundArea);
print("");

bPhrase2 = 1;
foundArea = integrate(minX,maxX,0);
print("bPhrase2: ",bPhrase2);
print(foundArea);
print("");

bPhrase2 = 10;
foundArea = integrate(minX,maxX,0);
print("bPhrase2: ",bPhrase2);
print(foundArea);
print("");

bPhrase2 = -1;
foundArea = integrate(minX,maxX,0);
print("bPhrase2: ",bPhrase2);
print(foundArea);
print("");

bPhrase2 = -10;
foundArea = integrate(minX,maxX,0);
print("bPhrase2: ",bPhrase2);
print(foundArea);
print("");
*/














